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We determine the timescales associated with turbulent diffusion and isotropization in closure models using anisotropically 
forced and freely decaying turbulence simulations and to study the applicability of these models. We compare the results 
from anisotropically forced three-dimensional numerical simulations with the predictions of the closure models and obtain 
the turbulent timescales mentioned above as functions of the Reynolds number. In a second set of simulations, turning the 
forcing off enables us to study the validity of the closures in freely decaying turbulence. Both types of experiments 
suggest that the timescale of turbulent diffusion converges to a constant value at higher Reynolds numbers. Furthermore, 
the relative importance of isotropization is found to be about 2.5 times larger at higher Reynolds numbers than in the more 
viscous regime. 



^7 1 Introduction 



The dynamics of many astrophysical large-scale flows such 
as solar and stellar differential rotation are strongly con- 
trolled by velocity correlations at smaller scales. These cor- 
relations are referred to as components of the Reynolds 
stress tensor. It is well known that in rotating stratified con- 
vection the Reynolds stress tensor is anisotropic (Kippen- 
hahn 1963), which then leads to the generation of differen- 
tial rotation (Rudigeret al. 1980, 1989). The Reynolds stress 
is defined as the average of products of components of ve- 



locity fluctuations, i.e., Ri 



where u = U — U 



l 3 — % 3 ' 

is the fluctuation of the velocity U about its mean U . Here 
and in the following, overbars denote mean quantities, and 
for the purpose of this paper we shall restrict ourselves to 
volume averages. 

Of particular interest are the equations governing the 
evolution of iZy. In the astrophysical context, such model 
equations have been derived by Ogilvie (2003) and Garaud 
& Ogilvie (2005); see also Kapyla & Brandenburg (2008), 
Snellman et al. (2009), and Garaud et al. (2010). Such equa- 
tions contain all the linear effects such as shear and rota- 
tion exactly. They usually also contain a driving term, Fij, 
through which energy is injected into the system, as well 
as viscous and turbulent damping terms. Finally, there often 
is a term that describes, in a somewhat more ad-hoc fash- 
ion, the return to isotropy (Rotta 1951). The latter is impor- 
tant if the off-diagonal components happen to be different 
from zero due to some statistical perturbation. At least at 
the level of a thought experiment, one might ask how the 
system returns to isotropy after the effects that produced the 
anisotropy, e.g., rotation and stratification via the A-effect, 



© 201 1 WILEY-VCH Verliig GmbH & Co. KGaA, Wcinhcim 

have been turned off. Mathematically, the turbulent damp- 
ing corresponds to terms involving triple correlations of the 
velocity while the term describing the return to isotropy 
comes from the interaction between components of veloc- 
ity and those of gradients of the pressure with the veloc- 
ity (Canuto 2009). Thus, in the absence of large-scale shear 
flows, rotation, gravity, or magnetic fields, we have 

Rij =Fij- r^Rij - T;^ 1 (Rij -±5ijR), (1) 

where the dot denotes a time derivative, R = Ra is the 
trace of R^, while r and Ti so are the relevant time scales 
describing turbulent diffusion and the return to isotropy. 

Two very similar ways to characterizing these 
timescales have been proposed, both of which as- 
sume proportionality to the eddy turnover time, and 
fcf is the wavenumber of the energy-carrying eddies. 
To = (urmsfcf)" 1 , where u rms is the rms velocity. In 
the standard minimal r-approximation (hereafter MTA) 
(Blackman & Field, 2002, 2003) the return to isotropy is 
not accounted for, and r is assumed constant in time. The 
value of r can be expressed in terms of To by defining a 
Strouhal number, St, via 



St ro. 



(2) 
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If the isotropization term is included in MTA, t- 1so is, like 
r, also considered constant. In an approach used by Ogilvie 
(2003), the rms velocity is written as u ima = i? 1//2 , and 
dimensionless fit parameters are introduced to quantify r 

and T iso : 

r- 1 = dhR 1 ' 2 , rrj = c 2 kiR 1 ' 2 . (3) 

Besides the non-vanishing isotropization term, the main dif- 
ference between these models is the nature of the eddy 
turnover time: in MTA it is usually constant, while in the 
Ogilvie approach it depends on the local and instantaneous 
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value of R. The latter model can be thought of as an exten- 
sion of the former to the case where u lms varies. 

There seems to be some diversity regarding the recom- 
mended choice of the coefficients c\ and 02- For the ra- 
tio C1/C2, Garaud & Ogilvie (2005) found the value 0.67, 
while in the additional presence of magnetic fields, Ogilvie 
(2003) found 0.87, and Liljestrom et al. (2009) found 0.86. 
The work mentioned above has attempted to compute these 
coefficients as fit parameters in models where additional ef- 
fects such as shear, rotation, and gravity are present. Such 
effects may however distort the results for c\ and c%, which 
characterize effects that are present even without the afore- 
mentioned processes. 

A goal of this paper is to determine the two non- 
dimensional coefficients ci and C2 using direct numerical 
simulations (DNS). We compute c\ and C2 here by imposing 
an anisotropic forcing term such that certain off-diagonal 
terms of its correlation matrix are non- vanishing. We use 
two independent methods to estimate the parameters: firstly, 
by comparing the steady state values for R and Rij to the 
strength of the forcing, and secondly by observing the be- 
havior of the system once the forcing is turned off, that is 
freely decaying turbulence. The predictions of the MTA and 
the Ogilvie approach regarding the behavior of the system 
in the latter case differ from one other, thus allowing us to 
assess the assumptions behind the two closures. 

2 The model 

We consider here a fully compressible gas with an isother- 
mal equation of state for which the pressure p is propor- 
tional to the density p with p = pc^, where c s = const 
is the isothermal sound speed. The computational domain is 
assumed Cartesian x = (x,y, z) with triply periodic bound- 
ary conditions. In some of our decay calculations, we start 
from a run where the Coriolis force is included, which is 
characterized by the angular velocity vector f2 = (0, 0, il). 
The equation of motion and the continuity equation can then 
be written as 
DC/ 
~Dt 



1 



= -<Vlnp-2ftx U + f + -V ■ (2vpS), (4) 



Dlnp 
Dt 



(5) 



where D/Dt = d/dt + U ■ V is the advective deriva- 
tive, Sy = \{Uij + Uj t i) — ^Sij"V ■ U is the traceless 
rate of strain matrix, commas denote partial differentiation, 
t is the time, and v is the kinematic viscosity. The forcing 
term is an adaptation of a previously used (Brandenburg 
2001) isotropic nonhelical forcing expression, / lso , which 
is monochromatic with wavenumber fe, whose modulus lies 
in a narrow band around an average wavenumber fcf, and 
the forcing is <5-correlated in time such that k[(t) changes 
abruptly from one time step to the next. The isotropic forc- 
ing function is written as / = Nfke lk ^' x , where N is a 
normalization factor, and fk = exk (with random unit 



vector e) to ensure that the forcing is solenoidal. Both e 
and k are random and non-parallel to each other. Next, we 
introduce a finite xy correlation by writing the forcing term 



as 



(6) 



= f +v(xf y +yf x ), 

where x and y are unit vectors in the x and y directions, re- 
spectively, and a is a non-dimensional parameter measuring 
the degree of anisotropy. Note that 



My = (i + a 2 )frfy° + ^[(/r) 2 + {fr? 



(7) 



and since f x so f y so vanishes on the average, f x f y has a pos- 
itive definite mean. This then implies that in the Reynolds 
equations (fl]i the forcing tensor 

Fij = p{uifj + Ujfi) (8) 

is also anisotropic with F xy on the average. 

To compute the effective timescales we consider steady 
state conditions in which case Eq. (Q~|) implies 

r- 1 - (F)/(R), (9) 

with F = Fa being the trace of F^ , and 

T-'+T.-^iF^/iRxy), (10) 

where angle brackets now denote time averages. 

A relevant control parameter is the Reynolds number, 
defined as 

Re=^, (11) 



vkr 



which is varied between 3 and 200. In some of the decay cal- 
culations that are initialized with rotation, we used a Cori- 
olis number, Co = 2Cl/u Ims k{ of order unity. In all other 
cases we have Co = 0. 

3 Results 

We have produced three-dimensional DNS models with 
anisotropic forcing varying both the Reynolds number and 
also the effective wavenumber of the forcing, fcf. Firstly, we 
determine t _1 and r^ 1 by comparing the steady state val- 
ues for R and R^ to the strength of the forcing in Sect. 13.11 
In these experiments the numerical resolution is 256 3 mesh- 
points. Secondly, we determine the inverse relaxation time 
scales from freely decaying turbulence in Sect. 13.21 Here, 
the numerical resolution is 128 3 meshpoints. 

3.1 Anisotropically forced turbulence 

The inverse relaxation timescales r _1 and r^ 1 measured 
from anisotropically forced turbulence in a steady state with 
varying Reynolds number and effective forcing wavenum- 
ber are shown in Fig. Q] The results show a clear decline 
of t -1 and T-~^ toward larger values Re. At the same 
time, Tj~o is about 2.5 times larger than r _1 , implying that 
Cf/c2 ~ 0.4, which is somewhat smaller than the values 
quoted in the literature; see Sect. Q] 
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Fig. 1 Dependence of the inverse relaxation time scales 
(normalized by the dynamical value t^ 1 on Re. Solid and 
dashed lines are for r _1 and tZ£, respectively. 



3.2 Decaying turbulence 

In this section we deterime the value of the timescales r 
and Ti SO and obtain another estimate for these parameters by 
studying freely decaying turbulence. We also compare the 
validity of the assumptions behind MTA and the Ogilvie 
approach, since the closures predict decay behaviors that 
are different in the two cases. By letting the turbulence first 
achieve a saturated state and then turning off the forcing in 
our DNS we get a time series that can be compared with 
the predictions of the closure models. From Eq. ([TJ we can 
easily derive the time evolution equation for the trace of 
the Reynolds tensor by summing over the diagonal com- 
ponents: 

R = F-t- 1 R 1 (12) 

where the summation causes the contribution from the 
isotropization term to vanish. Let the forcing be set to zero 
at t = to and let R(to) = R(°\ If t _1 is assumed con- 
stant in MTA, this approach predicts exponential decay. By 
integrating Eq. ( fT2b in this case we have 

R=R (0) e -(t-t )/T_ (13) 

The Ogilvie approach, however, predicts inverse square- 
type decay: 

-2 



R = 



1 



-cik f (t 



to) 



(14) 



By plotting Eqs. (13[ and (TBI) with the time series from 
DNS the behavior of the closures can be tested and the 
model parameters c\ and r estimated. We have performed 
two sets of runs, the results of which are summarized in 
Table Q] In Set F, we use the forcing scheme described in 
Sect. |2 while the runs in Set L were made using anisotropic, 
nonhelical forcing in combination with rotation (O ^ to) 
produce off-diagonal Reynolds stress components through 
the A-effect; see Kapyla & Brandenburg (2008) for a de- 
tailed description). The values listed in the table were ob- 
tained by fitting Eqs. ( fl~3l l and (TBt to the DNS results. Two 



Table 1 The model parameters c\, c%, t" 1 and t-~^ ob- 
tained from the DNS of freely decaying turbulence. The su- 
perscripts b and I refer to the beginning and the late parts of 
the time series. 



Run kf/ki 


Re 






A 


— 1 / —1 

r iso /To 


C b 2 


4 


LI 


3 


53 


0.11 


0.12 


0.16 


0.08 


0.14 


_ 


L2 


3 


55 


0.11 


0.12 


0.16 


_ 


_ 


_ 


L3 


3 


61 


0.12 


0.14 


0.16 


_ 


_ 


_ 


L4 


3 


79 


0.12 


0.14 


0.17 


0.04 


0.07 


_ 


L5 


1.5 


147 


0.08 


0.09 


0.15 


_ 


_ 




L6 


10 


14 


0.19 


0.23 


0.27 


0.06 


0.08 




L7 


10 


32 


0.13 


0.15 


0.18 








L8 


3 


113 


0.11 


0.13 


0.16 


0.08 


0.08 




L9 


3 


191 


0.10 


0.11 


0.15 


0.06 


0.10 




Fl 


3 


24 


0.19 


0.21 


0.25 


0.13 


0.15 




F2 


3 


53 


0.13 


0.13 


0.19 


0.19 


0.23 


0.05 


F3 


3 


92 


0.13 


0.14 


0.18 


0.27 


0.27 


0.07 


F4 


1.5 


55 


0.13 


0.15 


0.25 


0.32 


0.32 




F5 


1.5 


115 


0.12 


0.13 


0.22 


0.31 


0.31 


0.12 


F6 


1.5 


192 


0.12 


0.13 


0.20 


0.10 


0.11 


0.04 


F7 


10 


5 


0.44 


0.48 


0.65 


0.09 


0.23 




F8 


10 


13 


0.23 


0.27 


0.32 


0.13 


0.15 




F9 


10 


24 


0.16 


0.18 


0.24 


0.17 


0.17 





examples of such a fit can be seen in Fig. [2] The solid lines 
represent the DNS data, the dashed red lines the decay be- 
havior predicted by the MTA. The yellow and blue dotted 
lines are the corresponding prediction of the Ogilvie closure 
with two different values for c\, denoted with c\ and c\ for 
the determination of which the beginning and later parts of 
the DNS time series was used, respectively. The two alter- 
native fits for the latter model have been introduced because 
of the changing nature of the process. As we can see, the 
decay generally follows the exponential pattern at first, but 
in the later stages power-law behavior similar to the predic- 
tion of the Ogilvie model takes place. However, eventually 
the DNS results move away from both predictions. 

This kind of changing behavior is observed in all of the 
decay models, and the temporal span of the validity of var- 
ious predictions vary between the runs. This can be seen in 
Fig. |2 in which the upper panel shows the fit to the DNS 
data from Run F7, and the lower panel shows a correspond- 
ing fit to the data from Run F9: while the exponential pre- 
diction of MTA seems to apply for approximately the same 
duration in both panels, the Ogilvie approach has clearly 
a different range of applicability. Table Q] lists the different 
fit parameters c\, c[ and r _1 /t^ 1 obtained from the decay 
models. The values for c\ are generally very close to the 
values of r _1 /tq , while c\ tend to be somewhat larger. 
Actually, if one puts c\ = r _1 /r^ 1 , the resulting curve has 
MTA prediction as a tangent at to- 

Parameters t !SO and c-2 can be estimated by studying 
the decay of the off-diagonal components of the Reynolds 
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255 260 265 270 
t/r 




360 380 400 420 

t/r 



Fig. 2 The time evolution of R in freely decaying turbu- 
lence. Dotted and dashed lines show the decay predictions 
of the Ogilvie model and MTA, respectively, with suitable 
values for c\ and r , and the solid line is the DNS time se- 
ries. The upper panel shows the results from Run F7 and 
lower panel results from Run F9. 




360 380 400 420 440 460 480 500 
t/r 



Fig. 3 The time evolution of R xy in freely decaying tur- 
bulence. Dotted and dashed lines show the decay prediction 
of the Ogilvie model and MTA, respectively, with suitable 
values for c 2 and t , and the solid line is the DNS time se- 
ries. The runs displayed are the same as in Fig. [2] 



stress. The time evolution equation for Rij in the forced 
non-diagonal case reads 



Rij Fij [T 



)Ri 



(15) 



Now, let R^ (to) — R\y . Assuming Tj SO constant in the case 



of MTA we have again exponential decay: 

R ij =RS ) e- (t - toKT " +T ^ ) - 



(16) 

To get the corresponding result for the Ogilvie model one 
needs to use Eq. ( TBI to solve for ^/~R and integrate over 
time. The final result reads 

o c l+ c 2 

1 + — =— cifcf(*-*o) • (17) 



Rij - R %j 



The DNS results are compared with the predictions from 
the closure models in Fig. [3] Again we show two alternative 
versions for the behavior of the Ogilvie model with differ- 
ent values for c 2 , c 2 and c 2 , with the same reasoning as with 
c\. According to Eqs. ([Tol l and ( fT7b , the decay of R xy de- 
pends on the relaxation parameters t and c\ as well as the 



dedicated isotropization parameters Tj SO and c 2 . Using the 
estimates for the relaxation terms obtained from the decay 
of R we can determine the isotropization terms by treating 
them as the only free parameters of the models and finding 
a reasonable fit, like before. In the case of c 2 we have used 
the initial value c\ for this purpose. 

The results for the isotropization terms are summarized 
in TableQ] A problem in many runs is that the fluctuations of 
the off-diagonal components of the Reynolds stresses can be 
larger than their average value, causing their sign to change 
frequently. In the decay phase the time series of these runs 
tend to contain strong oscillations right from the beginning. 
The oscillations are similar to what can be seen in Fig. [3] 
and they make finding an unambiguous fit very challeng- 
ing. In some cases a suitable fit would have required nega- 
tive values for the parameter c 2 . For these cases, no value 
is given in Table Q] This problem manifests itself mostly 
in Set L. Thus, the most reliable results come from Set F, 
where R xy get non-zero mean values more consistently, and 
fluctuations are not too large. We see that T;7 /tq 1 an d c \ 
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Fig. 4 Inverse relaxation time scales (normalized by the 
dynamical value Tq 1 ) as functions of Re obtained from 
the decay models. Dotted and dashed lines are for t _1 
and Tj~; respectively. The diamonds represent runs with 
kt/ki = 1.5, triangles fcf/fei = 3 and asterisks fcf/fei = 10. 



obtain very similar values, while c l 2 is mostly very small or 
zero. Equation (TTTb implies that with c 2 = the decay of the 
off-diagonal components of the Reynolds stresses should 
behave like the decay of R described by Eq. (fT4] i. so the 
vanishing of c\ may indicate the isotropization switching 
off. But then again, it is seen in Fig. [3] that even with van- 
ishing C2 the prediction becomes gradually worse as time 
progresses, and in the lower panel the period of validity is 
restricted to a brief intersection. Large fluctuations are an- 
other source of ambiguity near the end of the time series. 

Figure|4]contains the same results as Fig.Q] but obtained 
for the decay models. Due to the ambiguity of the results 
from the Set L, only results from Set F are shown for Tj~ Q . 
In both figures the overall trend is similar: r _1 is large with 
small Reynolds numbers, and decreases as Re increases. 
Unlike in Fig.Q] in Fig.^r^ 1 generally increases with in- 
creasing Re, and eventually becomes greater than r _1 . It 
would seem that the results for r _1 approach some con- 
stant value at high Reynolds numbers, but more simulations 
with higher Reynolds numbers would be needed to verify 
this. Increasing t -1 with decreasing Re may explain, why 
the nature of the decline of R changes in the decay models. 
If we take u rms = i? 1 / 2 in the decay phase, the effective 
Reynolds numbers falls accordingly. This would mean that 
r -1 changes during the simulation, leading to a different 
behavior. 



different closure model predictions, namely the minimal tau 
approximation and the Ogilvie approach. 

Our results from the steady-state forced turbulence mod- 
els show that the values of r _1 , describing the diffusion pro- 
cess, and Tj~; , , describing the isotropization process, depend 
on Re for small and intermediate values, but show clear 
signs of convergence for larger values. In particular, it turns 
out that Tj7 is clearly larger than t _1 , and that their inverse 
ratio is around 0.4, which is somewhat less than the results 
published earlier in the literature. 

Our models of freely decaying turbulence show that, 
while the decay is exponential at first, as predicted by the 
MTA with a constant r, it deviates from this pattern in 
the later stages, following a power-law behavior much like 
the one predicted by the Ogilvie approach. Finally also the 
Ogilvie prediction breaks down far away from the switch- 
off point of the forcing. 
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4 Conclusions 

In this study we have investigated anisotropically forced hy- 
drodynamic turbulence, and determined the timescales re- 
lated to the diffusion and isotropization processes from our 
DNS models. The obtained results were compared to two 
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